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1  Objectives 

The  space  situational  awareness  (SSA)  mission  encompasses  intelligence,  surveillance  of  all  space  objects, 
and  the  prediction  of  space  events,  threats,  and  activities.  Specifically,  the  mission  requires  knowledge  of 
the  object  trajectory  (orbit)  and  type  (active  satellite,  decommissioned  satellite,  debris,  etc.)  for  all  objects  in 
the  near-Earth  environment.  For  objects  that  arc  currently  under  active  control,  one  would  also  like  to  know 
their  current  activities,  capabilities,  and  expected  future  actions.  These  types  of  higher-level  knowledge 
must  be  supported  by  a  significant  amount  of  raw  data  collection.  Space  surveillance  is  that  component  of 
SSA  focused  on  the  detection  of  space  objects  and  the  use  of  multisource  data  to  track  and  identify  space 
objects.  Currently,  it  can  take  weeks  to  establish  correct  orbits  on  objects  of  interest.  Our  approach  to  this 
problem  is  to  use  a  full  multiple  hypothesis  tracking  (MHT)  correlation  algorithm  developed  under  previous 
AFOSR  grants,  but  modified  to  accommodate  space  surveillance.  This  new  algorithm  will  be  called  multiple 
hypothesis  correlation  (MHC). 

While  manual  and  rudimentary  automated  correlation  schemes  may  suffice  for  known  orbits  or  for  a 
small  number  of  uncorrelated  tracks  (UCTs),  an  automated  system  will  be  required  in  the  future  to  handle 
an  anticipated  substantial  increase  in  the  number  of  space  objects  arising  from  space  events  and  higher- 
resolution  sensors.  Core  challenges  in  the  development  of  such  an  automated  system  include  dense  target 
environments,  the  need  to  establish  robust  tracks  in  a  timely  fashion,  limited  data  with  large  coverage  gaps, 
system  biases  and  residual  biases  and  their  impact  on  orbit  state  propagation,  and  unresolved  closely  spaced 
objects. 

From  a  multiple  target  tracking  perspective,  our  objectives  in  this  program  arc  to  initiate  orbit  trajectories 
on  newly  observed  objects  as  quickly  as  possible;  to  improve  our  ability  to  maintain  continuous  tracks  in 
the  presence  of  large  coverage  gaps;  and  to  piece  together  UCTs  from  different  sensors.  To  meet  these 
objectives,  research  is  required  in  the  following  three  algorithm  components. 

1.  Uncertainty  Modeling  and  Management.  A  prerequisite  to  a  statistically  rigorous  and  information 
theory-based  approach  to  tracking  objects  in  space,  sensor  resource  management,  and  conjunction 
analysis  is  the  consistent  characterization  of  the  uncertainty  exhibited  by  a  stochastic  state.  In  space 
surveillance,  the  challenges  to  achieving  this  objective  arc  the  nonlinear  dynamics,  long-term  prop¬ 
agations,  and  spai'sc  data  environment  requiring  the  development  of  advanced  methods  such  as  the 
adaptive  Gaussian  sum  filter  and  the  sliding  window  batch  estimation  filter.  The  status  of  the  current 
research  towards  this  objective  is  detailed  in  Subsection  1.1.1. 

2.  Multiple  Hypothesis  Tracking.  The  data  association  or  correlation  problem  is  that  of  determining 
which  reports  (tracks,  tracklets,  measurements,  features)  emanate  from  which  object.  For  widely 
spaced  objects,  reports  are  assigned  to  objects  if  they  arc  within  some  gate  (uncertainty  region)  of 
the  object.  (The  method  is  called  nearest  neighbor.)  On  the  other  hand,  multisensor  association 
(correlation)  in  a  dense  target  environment  is  a  challenging  problem  and  multiple  hypothesis  tracking 
(MHT)  techniques  offer  the  best  potential  solution.  For  marginally  detectable  or  dim  targets,  MHT 
can  be  augmented  by  track-before-detect  algorithms  that  output  detections,  but  not  formal  tracks,  or 
possibly  enhanced  resolution  techniques.  Thus,  what  is  needed  is  an  MHT  that  can  adapt  to  the 
complexity  of  the  problem  thereby  covering  the  range  between  widely  spaced  and  closely  spaced 
objects  as  well  as  luminous  and  dim  targets.  Perhaps  one  of  the  most  overlooked  problems  in  the  use 
of  multisource  data  and  modern  data  fusion  and  correlation  algorithms  is  that  of  sensor  biases.  A 
major  new  necessity  is  joint  bias  estimation  and  correlation,  which  may  be  needed  for  overlapping 
coverage,  but  is  mandatory  in  handover  situations,  especially  in  large  coverage  gaps.  Just  as  it  is 
essential  to  capture  the  uncertainty  in  state  estimation  via  nonlinear  filtering  techniques,  it  is  also 
essential  to  capture  the  uncertainty  in  the  association  process.  We  call  this  association  ambiguity 
and  we  achieve  this  objective  by  producing  the  “probability  of  association.”  Such  an  assessment  is 
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particularly  important  for  processing  UCTs  and  closely  spaced  objects.  The  status  of  the  current 
research  tow  ards  this  objective  is  detailed  in  Subsection  1.1.2. 

3.  Unresolved  Closely  Spaced  Objects.  Sensor  resolution  is  a  major  obstacle  in  establishing  orbits 
on  individual  objects  with  limited  data  in  a  timely  fashion;  thus,  it  becomes  necessary  to  enhance 
the  resolution  of  sensors  and  extract  consistent  sets  of  measurements  from  which  individual  orbits 
on  individual  objects  may  be  established.  The  status  of  the  current  research  tow  ards  this  objective  is 
detailed  in  Subsection  1.1.3. 

While  the  above  objectives  remain  a  focus,  Numerica  continues  the  work  on  network-centric  tracking 
that  was  initially  started  with  AFOSR  funding  and  has  now  transitioned  to  a  Phase  II  SBIR  at  AFRL/SNAT 
and  a  Phase  II  SBIR  with  Department  of  the  Army  with  a  transition  path  to  the  SIAP  Joint  Program  Office. 

1.1  Status  of  Effort 

The  status  of  the  current  research  program  is  summarized  in  this  section.  The  individual  objectives  outlined 
above  arc  each  addressed  in  the  following  subsections. 

1.1.1  Uncertainty  Modeling  and  Management 

In  this  subsection,  we  review  the  general  Bayesian  framework  for  nonlinear  filtering  and  describe  the  most 
promising  methods  for  achieving  correct  uncertainty  consistency  including  the  Gaussian  sum  filter  and  the 
sliding  window  batch  estimation  filter. 

1.1. 1.1  Nonlinear  Estimation  and  Filtering  The  problem  of  nonlinear  filtering  requires  the  definition 
of  dynamical  and  measurement  models.  It  is  assumed  that  the  dynamic  state  x(t)  6  Mn  at  time  t  evolves 
according  to  the  continuous-time  stochastic  model, 

x'(t)  =  f(x(t),t)  +  G  (x(t),t)w(t),  (1) 

where  /  :  Mn  xR^  Mn,  G  :  Mn  x  K  — ■>  Mnxm,  and  w(t)  is  an  m-dimensional  Gaussian  white  noise 
process  having  spectral  density  matrix  Q (t)  with  E[w(f)  w(t)t]  =  Q (t)  S(t  —  r),  t  >  r.  In  particular-, 
in  (1),  the  function  f  encodes  the  deterministic  force  components  of  the  dynamics  (e.g.,  gravity,  drag,  etc.) 
while  the  process  noise  term  models  the  stochastic  acceleration.  In  many  tracking  applications,  it  is  often 
convenient  to  work  with  a  discrete-time  formulation  of  the  dynamical  model  which  assumes  the  form 

Xk+l  f  T  Gr k{Xk)  UTfc,  (2) 

where  xk  =  x(tk),  fk  :  Rn  — ►  Mn,  G^  :  Mn  — >  Mnxm,  and  { w k}  is  an  m-dimensional  zero-mean  Gaussian 
white  noise  sequence  possessing  a  spectral  density  matrix  Q/.  e  Mmxm  such  that  E[io/c  wj]  =  Q/,  5k j-  In 
space  surveillance,  the  process  noise  term  is  often  very  small  and  discarded.  In  such  situations,  the  function 
fk  is  just  the  solution  flow  corresponding  to  the  continuous  model  (1)  with  w(t )  =  0. 

A  sequence  of  measurements  Zk  =  {z i, . . . ,  zk}  is  related  to  the  corresponding  kinematic  states  xk 
via  measurement  functions  hk  :  Mn  — >  M1'  according  to  the  discrete-time  measurement  model 


zk  =  hk(xk)  +  vk.  (3) 

In  this  equation,  {vk}  is  a  p-dimensional  zero-mean  Gaussian  white-noise  sequence  with  E[vkvJ]  = 
R k5kj.  More  general  filter  models  can  be  formulated  from  measurement  models  with  non-Gaussian  or 
correlated  (e.g.,  colored)  noise  terms  and  sensor  biases. 
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Figure  1:  The  predictor-corrector  step  for  the  recursive  Bayesian  state  estimator. 


In  the  Bayesian  approach  to  dynamic  state  estimation  [1],  one  constructs  the  posterior  probability  den¬ 
sity  function  (PDF)  of  the  state  based  on  information  of  a  prior  state  and  received  measurements.  Encap¬ 
sulating  all  available  statistical  information,  the  posterior  PDF  p(xk\Z^)  may  be  regarded  as  the  complete 
solution  to  the  estimation  problem  and  various  optimal  state  estimates  can  be  computed  from  it. 

In  a  recursive  filtering  approach,  measurement  data  is  processed  sequentially,  rather  than  as  a  batch. 
Given  the  initial  density  of  the  state  p(x o)  =  p{xq\Zq),  the  PDF  p(xk\Z^ )  is  obtained  recursively  in  two 
stages,  namely  prediction  and  correction,  as  illustrated  in  the  flowchart  of  Figure  1.  In  the  case  of  discrete¬ 
time  dynamics,  the  former  is  obtained  from  the  transitional  density  p(xk |aJfc_i)  in  conjunction  with  the 
Chapman-Kolmogorov  equation  defined  in  the  prediction  step  box  of  Figure  1 .  In  the  case  of  continuous¬ 
time  dynamics,  the  time  evolution  of  the  predicted  PDF  p  =  p(x,  t]Zk~i)  for  t  >  tk-i  is  governed  by  the 
Folclcer-Planck-Kolmogorov  equation  (FPKE)  [2] 

%  =  +  ^tr[V*v£(pGQGT)],  (4) 

where  \7X  is  the  gradient  with  respect  to  x  viewed  as  a  column  operator.  In  the  measurement  update  stage, 
also  called  the  correction  or  fusion  step,  the  density  p{zj;\x},)  is  evaluated  from  the  measurement  model  (3). 
The  term  p(zji\ Z^-i)  in  the  denominator  of  the  correction  step  is  called  the  prediction  error  and  appears  in 
the  likelihood  ratios  for  scoring  an  assignment  of  a  report  (i.e.,  z/J  to  a  track  (see,  for  example,  Poore  [3]). 
Thus,  its  accurate  evaluation  is  critical  for  correct  data/track  association  (correlation). 

Analytical  solutions  to  the  FPKE  and  to  the  prediction  and  correction  equations  in  Figure  1  are  generally 
intractable  and  are  only  known  in  a  few  restrictive  cases.  In  practice,  models  are  nonlinear  and  states  can 
be  non-Gaussian;  one  must  be  content  with  an  approximate  or  suboptimal  algorithm  for  the  Bayesian  state 
estimator.  While  the  extended  Kalman  filter  (EKF)  and  unscented  Kalman  filter  (UKF)  are  used  extensively 
in  ah'  and  missile  tracking,  they  only  represent  state  uncertainties  by  a  covariance  matrix  and  this  need  not 
be  adequate  in  the  space  surveillance  environment.  Because  of  the  need  to  propagate  uncertainties  over  ex¬ 
tended  time  intervals  in  the  absence  of  measurement  updates,  higher-order  cumulants  (e.g.,  skewness,  excess 
kurtosis)  can  become  non-negligible  and  must  be  accounted  for  in  order  to  achieve  uncertainty  consistency. 

In  this  AFOSR  program  and  in  another  STTR  funded  by  AFOSR  (topic  AF09-BT1 1),  we  developed  a 
variety  of  nonlinear  filters  specialized  from  the  general  Bayesian  state  estimator  and  showed  that  many  of 
them  correctly  capture  statistics  beyond  a  Gaussian  state  and  covariance  and  retain  uncertainty  consistency 
over  long  propagations.  Table  1  summarizes  some  of  the  advantages  and  disadvantages  of  these  filters. 
Although  no  special  attention  was  given  to  the  extended  Kalman  filter  (EKF)  and  particle  filter,  they  are 
nevertheless  included  in  the  table  for  completeness.  We  do  not  propose  their  use  in  space  surveillance 
applications.  The  UKF,  while  costing  about  the  same  to  run  as  the  EKF,  is  more  accurate  and  numerically 
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Table  1:  Comparison  of  some  nonlinear  filters. 

Filter 

Advantages 

Disadvantages 

Extended  Kalman  Filter 

Ubiquitous 

Uncertainty  is  only  represented  by  a  co- 
variance;  partial  derivatives  of  the  dy¬ 
namics  and  measurement  models  are  re¬ 
quired;  generally  less  accurate  than  the 
UKF 

Unscented  and  Gauss- 

Hermite  Filters 

Derivative-free 

Uncertainty  is  only  represented  by  a  co- 
variance 

Edgeworth  Filter 

Can  represent  any  uncertainty  within  a 
desired  accuracy  by  increasing  the  num¬ 
ber  of  cumulants  in  the  representation; 
derivative-free 

Computationally  feasible  only  for 
weakly  non-Gaussian  densities 

Reduced  Edgeworth  Filter 

As  in  the  Edgeworth  filter,  but  op¬ 
timized  to  account  for  possible  near 
Gaussianity  in  one  or  more  state  space 
coordinates 

As  in  the  Edgeworth  filter 

Adaptive  Gaussian  Sum 
Filter 

Can  represent  any  uncertainty  within 
a  desired  accuracy  by  increasing  the 
number  of  Gaussians  in  the  mixture; 
straightforward  to  implement  by  run¬ 
ning  a  bank  of  UKFs  in  parallel 

Representation  of  weakly  non-Gaussian 
densities  is  expensive;  computationally 
feasible  only  in  a  few  dimensions 

State  Transition  Tensor 

Filter 

Can  represent  any  uncertainty  within  a 
desired  accuracy  by  increasing  the  order 
of  the  Taylor  series  expansion;  very  ef¬ 
fective  for  representing  both  weakly  and 
highly  non-Gaussian  densities 

Partial  derivatives  of  the  dynamics  and 
measurement  models  are  required;  inte¬ 
gration  of  the  STT  ODEs  is  expensive 

Particle  Filter 

Can  represent  any  uncertainty  within  a 
desired  accuracy  by  increasing  the  num¬ 
ber  of  “particles”  in  the  representation 

Very  expensive;  accuracy  scales  as 
0{N~1^2)  where  N  is  the  number  of 
particles 

more  stable.  The  particle  filter  is  usually  reserved  as  a  “last  resort”  and  would  be  too  expensive  to  use  in  an 
operational  system. 

Figure  2  shows  how  the  root-mean-square  (RMS)  cost  error,  a  metric  for  quantifying  uncertainty  con¬ 
sistency,  scales  according  to  the  dimension  of  the  representation  (i.e.,  the  number  of  degrees  of  freedom 
required  to  store  the  state  PDF)  and  the  number  of  sigma  points  (i.e.,  the  number  of  orbital  propagations 
required  in  the  filter  prediction  step).  The  “ideal”  filter  would  achieve  a  negligible  RMS  cost  error  with 
minimal  computational  cost  either  in  terms  of  the  dimension  of  the  represented  PDF  or  the  number  of  sigma 
points.  The  adaptive  Gaussian  sum  filter  (AGSF)  and  the  state  transition  tensor  filter  (STTF)*  show  the  most 
promise  for  accurately  representing  a  state’s  uncertainty.  However,  beware  that  the  STTF  results  are  based 
only  on  unperturbed  Keplerian  dynamics  in  which  the  implementation  of  the  filter  is  trivial.  It  remains  to 
see  how  the  STTF  performs  under  a  full  nonlinear  gravity  model. 

*The  state  transition  tensor  filter  is  not  based  on  the  propagation  of  sigma  points,  so  it  does  not  appear  in  the  right  figure. 
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Figure  2:  Dependence  of  the  root-mean-square  (RMS)  cost  error  on  the  dimension  of  the  representation 
(left)  and  the  number  of  sigma  points  (right).  Blue  disks  (•)  imply  use  of  the  EGM96  degree  and  order  70 
gravity  model  while  red  squares  (■)  imply  use  of  unperturbed  Keplerian  dynamics. 


1. 1.1.2  Gaussian  Sum  Filters  A  Gaussian  sum  is  a  mixture  density  of  the  form 


N 

p(x)  =  ^2waAf(x-,tia,PQ), 

OL=  1 


where  the  weights  wa  are  non-negative  scalars  which  sum  to  unity  and  Af(x:  p.  P)  denotes  the  Gaussian 
PDF  with  mean  /x  and  covariance  P;  i.e., 


1 

Y /  det(27rP) 


exp 


(5) 


The  Gaussian  sum  filter  (GSF)  is  based  on  a  fundamental  result  of  Alspach  and  Sorenson  [4]  which  states 
that  any  PDF  can  be  approximated  arbitrarily  close  (in  the  L1  sense)  by  a  weighted  sum  (mixture)  of  Gaus¬ 
sian  PDFs  henceforth  called  a  Gaussian  sum.  Thus,  Gaussian  sums  provide  a  mechanism  for  modeling 
non-Gaussian  densities  and  for  more  accurately  approximating  the  solution  of  the  FPKE.  Computationally, 
the  GSF  has  the  added  advantage  of  being  parallelizable  since  filters  such  as  the  EKF  or  UKF  act  inde¬ 
pendently  on  each  component  Gaussian  in  the  prediction  and  correction  steps.  With  regards  to  updating  the 
weights  within  the  filter,  such  a  scheme  is  clearly  dictated  from  Bayes’  rule  following  a  measurement  or 
track  update  (fusion).  Flowever,  there  is  not  complete  agreement  on  how  the  weights  should  be  updated  (if 
at  all)  following  a  prediction. 

One  key  feature  of  the  new  GSF  developed  in  this  effort  and  communicated  in  the  paper  [5]  is  the  absence 
of  a  procedure  to  update  the  mixture  weights  following  the  implementation  of  the  filter  prediction  step. 
The  accuracy  and  consistency  of  the  propagated  uncertainty  is  ensured  by  representing  the  Gaussian  sum 
in  (traditional  or  alternate)  equinoctial  orbital  elements  with  component  means,  covariances,  and  weights 
initially  selected  (by  solving  an  L2  optimization  problem  offline)  such  that  the  (square -root  version  of  the) 
UKF  [6],  when  acting  in  parallel  on  each  component,  accurately  approximates  the  solution  of  the  FPKE.  The 
number  of  Gaussian  components  required  to  achieve  an  accurate  approximation  is  chosen  adaptively  based 
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Figure  3:  Depiction  of  a  single  Gaussian  and  its  Gaussian  sum  approximation  undergoing  a  nonlinear  trans¬ 
formation. 


on  the  length  of  the  propagation  time  and  the  initial  error  (standard  deviation)  along  the  radial  direction 
(semi-major  axis  coordinate).  Consequently,  by  representing  the  PDF  in  the  equinoctial  elements,  the  new 
algorithm  achieves  superior  computational  efficiency  because  it  only  requires  a  Gaussian  sum  along  one 
dimension. 

While  others  have  proposed  methods  for  adapting  the  weights  based  on  various  online  L 2  optimization 
criteria  [7-9],  we  have  found  that  the  application  of  these  methods  to  our  existing  GSF  does  not  improve  un¬ 
certainty  consistency,  but  rather  causes  the  accuracy  of  the  Gaussian  sum  representation  to  degrade  slightly 
[5].  We  attribute  these  findings  to  the  fact  that  our  GSF  is  already  very  efficient  because  it  solves  an  L2 
optimization  once  offline  which  is  used  thereafter  in  any  scenario  through  a  lookup  table.  In  what  follows, 
we  motivate  and  define  the  specific  optimization  problem  we  solve. 

The  underlying  optimization  problem  fundamental  to  our  GSF  requires  the  refinement  of  a  single  Gaus¬ 
sian  PDF  into  a  Gaussian  sum.  The  component  means,  covariances,  and  weights  of  the  mixture  are  chosen 
such  that  each  Gaussian  component  remains  Gaussian  (up  to  a  prescribed  accuracy)  when  propagated  by 
the  UKF  under  the  nonlinear  dynamics.  The  refinement  methodology  is  illustrated  in  Figure  3.  Under  a 
nonlinear  transformation,  a  Gaussian  (represented  by  the  thick  black  ellipse)  need  not  be  mapped  to  a  Gaus¬ 
sian  (e.g.,  the  level  surfaces  of  the  transformed  distribution  could  look  crescent-shaped).  However,  in  a 
sufficiently  small  neighborhood,  any  (smooth)  nonlinear  map  will  be  approximately  linear.  Consequently, 
Gaussians  with  smaller  covariances  (represented  by  the  colored  elliptic  disks)  remain  more  Gaussian  than 
those  with  larger  covariances  under  the  nonlinear  mapping.  Therefore,  a  Gaussian  sum  refined  by  approx¬ 
imating  each  constituent  Gaussian  by  a  finer  Gaussian  sum  will  exhibit  better  behavior  through  nonlinear 
transformations.  It  suffices  to  optimally  refine  the  unit  one-dimensional  Gaussian  J\T(x,  0, 1);  refinement  of 
a  multivariate  Gaussian  with  an  arbitrary  covariance  is  obtained  by  a  series  of  linear  transformations  detailed 
in  [5]. 

We  now  derive  a  Gaussian  sum  approximation  of  the  unit  one-dimensional  Gaussian  J\T  (x;  0, 1).  The 
problem  has  an  obvious  trivial  solution  which  perfectly  approximates  the  target,  namely  the  Gaussian  sum 
with  a  single  component.  However,  as  motivated  earlier,  the  idea  is  to  construct  a  Gaussian  sum  approxima¬ 
tion  whose  component  standard  deviations  <ja  are  smaller  than  some  fixed  value  o  <  1  so  that  the  Gaussian 
sum  more  accurately  represents  the  true  distribution  under  a  nonlinear  transformation.  This  requirement 
leads  to  a  constrained  nonlinear  optimization  problem.  That  said,  we  assume  the  target  Gaussian  sum  has 
the  form 

N 

Pgs(x)  =  ^waM{x\^aio2a), 

OL=l 

for  some  predetermined  length  N  >  1  and  define 


E  =  ^\\Af(x-  0, 1)  -  pgs(x)  \\2l2  =  X- 


N 


Af(x\0,l)  -  y^^waN{x]na,o2a) 


a=l 


dx. 


(6) 


The  weights  wa,  means  /x(1 ,  and  standard  deviations  oa  of  each  component  are  determined  from  the  solution 


UNCLASSIFIED 


6 


Numerica  Data.  Use  or  disclosure  of  data  on  this  page  is  subject  to  the  restrictions  on  the  cover  of  this  document. 


Number  of  Gaussians 


Figure  4:  Plots  of  the  Lv  error  between  the  unit  one-dimensional  Gaussian  and  its  Gaussian  sum  approxi¬ 
mation  obtained  using  the  suboptimal  optimization  algorithm. 


of  the  following  L2  optimization  problem: 


Minimize  E 


wi  cri,...,CT  N 

N 

Subject  to 

wa  =  1,  wa  >  0,  a  =  1, . 

a=l 

,.,N, 

(7) 

Pi  <  P2  <  •  •  •  <  Pat, 

Op  <  <7  <  1 7  a  =  1, . . . ,  N. 


Noting  the  identity 

J  n2,V2)&x  =  Nin-L -p2;0,Pi  +  P2), 

it  follows  that  (6)  reduces  to 

E  =  -w1  Mm  —  wTn  -\ - — ,  (8) 

2  4VF 

where 

(■ w)a  =  wa ,  (n)a  =  jV(pa;  0,  cr2  +  1),  (M)a/3  =  AA(pQ  -  p^;  0,  a2  +  cr|). 

The  minimization  of  (8)  over  the  individual  weights,  means,  and  standard  deviations  subject  to  the 
constraints  in  (7)  is  a  difficult  nonlinear  optimization  problem.  A  suboptimal  yet  computationally  tractable 
algorithm  which  approximates  the  solution  of  (7)  can  be  obtained  as  follows.  Specifically,  we  constrain 
each  of  the  component  Gaussians  to  have  a  common  standard  deviation  oa  =  o  <  1  and  propose  /kmg  the 
means  n„  according  to 

Pa  =  -m  +  a  (a  —  1),  (9) 

for  cr  =  1, ...  ,N,  where  m  >  0  and  N  =  [T  +  2m/cr].  Note  that  the  means  are  evenly  distributed  with 
the  left-most  and  right-most  components  located  at  x  ~  ±m.  In  order  to  ensure  adequate  accuracy  around 
the  tails,  we  propose  setting  m  =  4  if  o  >  otherwise  m  =  6.  With  these  additional  constraints,  the 
optimization  problem  (7)  reduces  to  a  quadratic  programming  problem  which  is  straightforward  to  solve 
using  elementary  methods. 

Figure  4  shows  plots  of  the  L1,  L2,  and  L°°  error  between  the  unit  one-dimensional  Gaussian  and 
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Figure  5:  Scatter  plots  of  the  component  weights  and  means  computed  using  the  suboptimal  algorithm  (i.e., 
minimizing  only  over  the  weights)  and  by  solving  the  full  (optimal)  optimization  problem  (i.e.,  minimizing 
over  the  weights,  means,  and  standard  deviations). 


its  Gaussian  sum  approximation  computed  using  the  suboptimal  approach  described  above.  Up  to  around 
N  ~  100,  the  three  L'p  errors  all  decrease,  but  they  plateau  to  around  10-8  for  N  100.  We  believe  this 
plateau  is  not  because  of  numerical  stability  issues  (the  matrix  M  in  (8)  has  a  condition  number  on  the  order 
of  105)  but  simply  because  we  at  e  fixing  the  standard  deviations  and  means  via  (9)  and  only  minimizing  over 
the  unknown  weights;  this  need  not  be  optimal.  Notwithstanding  these  comments,  the  simulation  studies 
considered  in  the  paper  [5]  use  this  (suboptimal)  refinement  scheme  with  N  as  high  as  1000  without  any 
adverse  effects. 

We  conclude  this  section  by  presenting  some  preliminary  solutions  of  the  full  optimization  problem 
(7).  We  have  found  that  the  numerical  conditioning  of  this  problem  worsens  as  the  target  length  N  of  the 
Gaussian  sum  increases  thereby  rendering  double  precision  floating-point  arithmetic  largely  inadequate.  To 
facilitate  an  accurate  solution,  we  have  used  the  Maple  ‘Optimization’  package  which  takes  advantage  of 
built-in  library  routines  provided  by  the  Numerical  Algorithms  Group  with  the  ability  to  call  them  within  an 
arbitrary-precision  software  floating-point  environment. 

Figure  5  shows  various  scatter  plots  of  the  component  weights  and  means  computed  using  two  differ¬ 
ent  methods.  The  first  method  uses  the  suboptimal  approach  which  fixes  the  means  according  to  (9)  and 
subsequently  optimizes  over  the  weights  only.  The  second  method  solves  the  full  optimization  problem  (7). 
Each  subplot  assumes  a  different  target  length  N  and  standard  deviation  a.  One  interesting  observation  is 
that  the  optimal  component  standard  deviations  lie  on  the  constraint  boundary;  i.e.,  aa  =  a  for  all  a.  Most 
importantly,  by  solving  the  full  optimization  problem,  the  locations  of  the  means  are  no  longer  uniformally 
distributed  and  the  resulting  L2  error  is  smaller  than  that  obtained  by  solving  the  suboptimal  problem.  In 
particular,  to  get  the  1 2  error  down  to  1CU '  requires  about  N  =  28  Gaussians  using  the  suboptimal  method 
but  only  about  N  =  9  Gaussians  when  solving  the  fid  I  optimization  problem!  We  emphasize  that  the  full 
problem,  although  very  expensive  to  solve,  can  nevertheless  be  done  offline.  One  can  generate  a  library  of 
Gaussian  sum  approximations  to  A f{x\  0, 1)  for  various  values  of  N  and  a. 

1.1. 1.3  Batch  Estimation  and  Sliding  Window  Batch  Estimation  Given  a  sequence  of  m  measure¬ 
ments  Zm  =  {zi, . . . ,  zrn\  at  times  t\, ... ,  trn  commensurate  with  the  model  (3),  the  objective  of  initial 
orbit  determination  (IOD)  and  batch  estimation  is  to  obtain  a  representation  of  the  (joint)  posterior  PDF 
p(x o, . . . ,  xm\Zm),  where  xj,.  denotes  the  dynamical  state  of  the  system  at  time  t/,.,  and  to  extract  mean¬ 
ingful  statistics  (e.g.,  mean,  covariance)  from  it  in  a  consistent  manner.  It  is  assumed  that  the  state  evolves 
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according  to  the  discrete-time  model 

Xk+i  =  fk(xk)  +  wk,  (10) 

where  { w k }  is  a  zero-mean  white  noise  sequence  with  E[wk  wj]  =  Qk  6kj.  The  following  independence 
assumptions  are  implied  between  the  prior  xo,  the  measurement  noise  sequence  {vk },  and  the  process  noise 
sequence  {wk}: 

E[*0  wj]  =  0,  E[*0  vk]  =  0,  E[vk  wj }  =  0. 

Appealing  to  Bayes’  rule  and  the  above  assumptions,  the  joint  posterior  PDF  is  derived  in  Jazwinski  [2, 
§5.3]  and  is  found  to  be 


m  m 

p(x  0  5  •  •  •  5  'Em  |  Z  m  )  =  cp0{x0)  Y[pWk{xk  -  fk-l(xk- 1))  n Pvk(Zk  -  hk(xk)),  (11) 

k=  1  k=  1 

where  c  is  a  normalizing  constant,  and  po  is  the  prior  PDF  of  the  state  x$  at  time  to-  Further,  in  (11),  the 
pWk  and  pVk ,  for  k  =  1, ...  ,m,  arc  the  respective  PDFs  of  the  process  and  measurement  noise  processes.  In 
practice,  they  are  often  assumed  to  be  Gaussian  with  zero  mean  and  covariances  of  Q/  and  R&,  respectively. 

The  posterior  PDF  (II)  is  the  complete  description  of  the  uncertainty  of  the  state  at  each  of  the  measure¬ 
ment  times.  In  practice,  a  finite  dimensional  representation  of  the  uncertainty  is  sought.  Thus,  the  emphasis 
of  the  batch  estimation  problem  is  on  how  statistical  information  can  be  extracted  from  (11)  consistently 
and  accurately.  Nonlinear  optimization  theory  provides  a  framework  for  computing  the  modal  trajectory  or 
maximum  a  posteriori  (MAP)  estimate  of  (11).  For  a  Gaussian  prior  with  xq  ~  N(x o,  Po)  and  Gaussian 
noise  processes  as  described  above,  the  modal  trajectory  is  obtained  by  solving  the  least  squares  or  batch 
problem 

(*o,  ■  ■  • ,  xm)MAP  =  arg  max  p(x0, ...,  xm\Zm ) 

1  1  m,  1  rn  (12) 

=  arg  min  -||*0  -  *o|||-i  +  x  V  \\xk  ~  fk- i(*fc-i)|lo-i  +  -  V  || zk  -  hk(xk) 

Xo,...,Xm  Z  0  Z  —  1  Z  x,k 

k= 1  k= 1 

In  initial  orbit  determination,  we  do  not  have  a  prior;  the  term  \\\xo  —  £Co||e-i  is  removed  from  the  above 
formulation. 

Methods  for  solving  nonlinear  least  squares  problems,  such  as  Gauss-Newton,  full  Newton,  and  quasi- 
Newton  updates,  along  with  globalization  methods  such  as  line  search  and  trust  region  methods  including 
Lcvenbcrg-Marq uardt  [10],  arc  efficient  and  mature  and  will  not  be  discussed  further  here.  In  the  astrody- 
namics  community,  such  solution  techniques  arc  called  differential  correction  methods.  In  any  nonlinear 
least  squares  problem  such  as  (12),  one  must  provide  the  solver  a  starting  guess  in  order  to  initiate  the  differ¬ 
ential  correction  method.  This  is  the  initial  orbit  determination  (IOD)  problem.  In  the  case  of  measurement 
data  from  a  single  radar  or  EO  sensor,  a  first  estimate  can  be  obtained  using  the  classical  methods  of  Lam¬ 
bert  or  Gauss  (see,  for  example,  [11]).  Additionally,  for  angle-only  observations,  a  recent  algorithm  due 
to  Gooding  [12]  has  shown  promise  for  IOD  scenarios  involving  both  ground-based  and  space-based  EO 
sensors  [13]. 

Seen  as  a  hybrid  between  batch  estimation  and  the  sequential  prediction-correction  filter  discussed  ear¬ 
lier,  the  sliding  window  batch  estimation  filter  (SWBEF)  processes  n  of  the  m  measurements  (through  a 
batch  process)  over  a  sliding  window.  The  SWBEF  is  the  full  batch  estimation  algorithm  when  the  sliding 
window  encloses  the  entire  measurement  sequence  (n  =  m)  and  is  the  traditional  predict-correct  sequential 
filter  when  n  =  1. 

The  use  of  a  SWBEF  can  be  very  effective  in  addressing  the  problems  of  anomaly  detection  and  UCT 
resolution  because  of  the  need  to  propagate  states  and  uncertainties  over  long  time  gaps  and  to  precisely 
evaluate  likelihood  ratios  used  in  the  association  problem.  In  particular,  in  the  problem  of  orbit  determi- 
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Figure  6:  Schematics  for  the  sequential  prediction-correction  filter,  sliding  window  batch  estimation  filter 
(of  length  n  =  1),  and  full  batch  estimation. 


nation  in  the  low  Earth  orbit  (LEO)  regime  (e.g.,  a  two  minute  track  generated  from  radar  measurements 
every  ten  seconds  apart),  it  is  observed  that  the  posterior  PDF  p{xk\Zk)  is  often  well  approximated  by  a 
single  Gaussian  in  equinoctial  orbital  elements  [14]  or  Cartesian  ECI  coordinates  [15].  However,  the  pre¬ 
dicted  PDF  p(xk+i\Zk.)  can  become  highly  non-Gaussian  if  the  propagation  time  is  long  or  the  uncertainty 
along  the  radial  direction  (semi-major  axis)  of  the  posterior  is  large.  One  can  expend  tremendous  compu¬ 
tational  resources  in  accurately  representing  the  predicted  PDF  (using,  for  example,  a  high  fidelity  GSF  or 
STTF).  Yet,  when  the  predicted  PDF  p{xk+\\Zk)  is  updated  with  a  new  report  z^+ 1,  the  resulting  density 
p(xk+i\Z k+i)  tends  to  collapse  to  something  which  is  nearly  Gaussian.  The  SWBEF  bypasses  the  expen¬ 
sive  representation  of  highly  non-Gaussian  densities  arising  from  the  prediction  step  and  instead  waits  until 
a  new  report  (measurement,  track)  becomes  available  and  then  performs  the  prediction  and  correction  step 
simultaneously.  Further  details  on  the  SWBEF  are  included  in  the  paper  [16]. 

1. 1.1.4  Status  A  number  of  papers  have  been  submitted  to  journals  and  conference  proceedings  on  un¬ 
certainty  modeling  and  management  including  ones  on  Gaussian  sums  [5,  17-19],  the  Edgeworth  filter  [19], 
the  sliding  window  batch  estimation  filter  [16],  and  the  batch  estimation  problem  [15].  Future  work  on 
this  objective  will  aim  to  mature  the  Gaussian  sum  filter  and  the  underlying  optimization  problem.  Other 
alternate  nonlinear  filtering  algorithms  will  also  be  pursued  including  the  state  transition  tensor  filter. 

1.1.2  Multiple  Hypothesis  Tracking 

Data  association  or  correlation  methods  for  multiple  target  tracking  divide  into  two  broad  classes:  single 
frame  methods  and  multiple  frame  methods.  The  single  frame  methods  include  nearest  neighbor,  global 
nearest  neighbor  based  on  a  two-dimensional  assignment  solver,  and  joint  probabilistic  data  association 
(JPDA)  [20].  For  many  widely  spaced  objects  and  a  clear  background,  nearest  neighbor  may  be  appropriate, 
especially  for  space  surveillance.  For  noisy  backgrounds,  JPDA  can  mitigate  misassociation  by  updating  the 
track  with  a  weighted  combination  of  all  the  measurements  within  a  gate  of  the  projected  track  state.  The 
most  successful  of  the  multiple  frame  methods  are  multiple  hypothesis  tracking  (MHT)  [21]  and  multiple 
frame  assignment  (MFA)  [3,  22].  MFA  is  an  optimization-based  formulation  of  MHT  and  is  now  considered 
to  be  the  industry  standard  replacement  of  MHT.  MHT/MFA  methods  mitigate  misassociation  by  the  ability 
to  hold  difficult  correlation  decisions  in  abeyance  until  additional  information  is  available,  as  well  as  provide 
an  opportunity  to  change  past  decisions  to  improve  current  decisions.  In  dense  tracking  environments  such 
as  breakups  or  geoclusters  with  closely  spaced  objects,  the  performance  improvements  of  multiple  frame 
methods  over  single  frame  methods  are  significant,  making  MHT/MFA  the  preferred  solution. 

The  challenges  of  space  surveillance,  such  as  large  coverage  gaps  and  both  sparse  and  dense  target 
scenarios,  require  investigation  of  major  extensions  to  current  MHT/MFA  architectures,  including  bias  esti¬ 
mation  and  ambiguity  assessment.  A  major  new  necessity  is  joint  bias  estimation  and  correlation,  which  may 
be  needed  for  overlapping  coverage,  but  is  mandatory  in  handover  situations,  especially  in  large  coverage 
gaps.  We  now  discuss  these  requirements  and  present  some  MHT  results  obtained  through  this  effort. 
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1. 1.2.1  A  Brief  Review  of  MHT/MFA  Multiple  hypothesis  tracking  normally  uses  a  moving  window 
over  a  set  of  data  frames.  The  concept  of  a  “data  frame”  is  central  to  most  modern  data  association  or 
correlation  algorithms.  For  initial  orbit  determination,  the  data  frames  correspond  to  sets  of  sensor  reports 
in  which  each  target  is  seen  at  most  once  in  each  frame.  For  track  maintenance  or  extension,  one  data  frame 
corresponds  to  a  set  of  system  tracks,  while  the  remaining  frames  arc  composed  of  sensor  reports. 

jj  j^/f3 

Let  Zi  =  { zf  }  d  k  denote  a  sequence^  of  noise-contaminated  measurements  of  length  Ml  from  sensor 

j  at  time  Note  that  as  the  sensor  index  j  and  time  index  k  are  varied  the  number  of  measurements  Mk 
can  vary  due  to  field  of  view  issues,  false  alarms,  missed  detections,  etc.-f  We  now  define  A  as  the  list  of 
time-sensor  index  pairs  (k,  j)  such  that  Zk  is  non-empty.  The  cardinality  of  A  is  N.  Let  ZN  =  {Z3k}_ 4 
denote  the  data  from  these  N  sensor  scans,  potentially  covering  multiple  times  as  well  as  multiple  sensors. 
The  central  data  association  problem  in  multitarget  and  multisensor  data  fusion  can  be  generally  posed  as 
[3] 

Prf  H\Zn) 

H*  =  arg  max  LH(N),  LH(N)  =  >  (13) 

where  H  denotes  a  partition  of  the  data  into  tracks  and  false  alarms,  H$  denotes  a  reference  partition  in 
which  all  reports  arc  declared  to  be  false  alarms.  7i  denotes  the  set  of  all  feasible  data  partitions  of  the  data 
ZN ,  and  H*  denotes  the  optimal  partition.  Typically,  one  replaces  the  N  frame  likelihood  ratio  in  (13)  with 
a  negative  log-likelihood  score 

ch{N)  =  -  In  Lh( N),  (14) 

which  changes  (13)  to  a  minimization  problem. 

1. 1.2.2  The  Multidimensional  Assignment  Problem  A  derivation  of  the  problem  and  the  equivalence 
with  MHT  using  an  iV-scan  back  approximation  has  been  previously  given  by  Poore  [3]  and  will  not  be 
repeated  here.  The  following  is  a  brief  summary  of  the  formulation,  but  not  the  derivation.  Consider  a 
layered  graph  Q  with  N  distinct  node  sets  (layers)  ,4/.  =  {0, 1, ... ,  my.}  indexed  by  keK  =  {i,...,N} 
and  arcs  A  C  A\  x  •  •  •  x  An.  (In  tracking  applications,  a  node  set  Ay.  is  called  a  frame  of  data.)  A  specific 
arc  in  A  is  denoted  by  the  iV-dimensional  multi-index  a  =  (a(l),  a(2), . . . ,  a(N))  where  each  a(k)  €  A^. 
Next,  define  a  section  Am  =  {a  6  A  \  a(k)  =  £}.  For  each  k  G  K,  associate  to  each  section  Am  a 
nonnegative  integer  hm  >  1  for  £  =  1, . . . ,  ?/;/,.  that  will  denote  the  number  of  times  the  index  (k.  £),  i.e., 
can  be  assigned. 

In  addition,  xa  is  assumed  to  be  a  zero-one  variable  for  each  a  G  ,4  with  at  least  two  nonzero  indices. 
The  iV-dimensional  problem  (multidimensional  assignment  problem  of  dimension  N)  can  then  be  expressed 
as 


Minimize  caxa 

a£A 

Subject  To  ^2  xa  <  nM  f°r  k  F  K  and  £  €  A^  \  {0},  (15) 

xa  G  {0, 1}, 

where  each  ri}.p  is  the  number  of  times  a  report  £  =  1, . . , ,  rrip.  on  frame  k  can  be  assigned.  In  most  tracking 
applications,  rikp  =  1;  however,  as  suggested  by  real  orbital  data  provided  by  JFCC  SPACE,  an  object  may 
be  seen  more  than  once  in  a  pass  over  a  sensor  and  thus  there  is  a  need  to  allow  multiassignment.  This  is 
a  new  feature  of  the  association  problem  for  space  surveillance.  The  cost  ca  is  the  aforementioned  negative 

1  No  distinction  is  made  between  objects  being  observed  as  the  association  (which  is  being  solved  for)  is  not  known  a  priori. 

4  In  space  tracking  false  alarms  are  rare,  but  missed  detections  are  quite  common. 
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log  of  a  likelihood  ratio  represented  above  as  the  ratio  of  two  probabilities  with  the  normalizing  probability 
that  each  track  let-  is  an  uncorrelated  track  (or  false  alarm  in  measurement  based  tracking). 

While  the  association  problem  can  be  formulated  as  a  zero-one  linear  programming  problem,  the  as¬ 
signment  formulation  has  been  shown  to  be  at  least  two  orders  of  magnitude  faster  on  larger  problem  sets. 
We  propose  the  use  of  four  classes  of  algorithms:  (a)  branch  and  bound  based  on  Lagrangian  relaxation  to 
produce  a  guaranteed  optimal,  (b)  a  Lagrangian  relaxation  algorithm  for  a  heuristic,  (c)  a  partial  branch  and 
bound  to  improve  the  Lagrangian  heuristic,  and  (d)  an  anytime  algorithm  that  can  control  memory  usage  and 
runtime.  For  (a)  and  (c)  we  make  use  of  A*-search  to  compute  the  k-best  solutions  for  the  multidimensional 
assignment  problem.  This  provides  the  basis  for  ambiguity  assessment. 

1. 1.2.3  Treatment  of  Biases  A  key  challenge  in  the  correlation  and  fusion  of  multiple  sensors  is  that  of 
dealing  effectively  with  sensor  biases  and  navigation  errors,  for  which  we  also  use  the  term  “biases.”  Ah 
fusion  systems  that  integrate  multiplatform  multisensor  data  require  that  sensor  registration  be  performed 
to  remove  biases  before  correlation  and  fusion.  Otherwise,  the  data  from  each  sensor  will  be  misaligned, 
and  performance  will  degrade.  While  “biases”  are  stochastic  in  nature,  they  represent  systematic  errors  that 
do  not  average  out  and  lead  to  misassociation  and  redundant  tracks  in  a  multisensor  environment.  Bias 
treatment  can  generally  be  divided  into  two  cases:  (i)  association  is  known  or  some  form  of  truth  is  available 
and  (ii)  association  is  unknown  and  no  truth  is  available. 

In  the  first  case,  the  issue  of  association  has  already  been  separated  out  and  the  bias  problem  can  be 
isolated.  This  is  the  case  when  the  association  between  reports  and  objects  is  known,  as  might  be  the  case 
with  widely  spaced  targets.  This  is  also  the  case  when  truth  (or  “fuzzy”  or  “approximate”  truth,  which  might 
include  some  small  noise)  is  available.  An  advantage  of  the  space  surveillance  tracking  environment  is  the 
availability  of  approximate  truth  data  for  a  small  collection  of  objects  with  well  known  trajectories.  These 
can  be  used  to  estimate  the  biases  for  the  sensors  in  the  network,  and  weekly  bias  averages  arc  often  used  to 
debias  the  following  week’s  data.  This  calibration  and  debiasing  of  the  sensor  reports  and  navigation  errors 
is  performed  using  nonlinear  least  squares  techniques,  which  can  also  reveal  the  observability  [23,  24].  The 
drawback  of  this  approach  is  that  the  sensor  biases  arc  not  constant,  and  weekly  updates  can  be  insufficient 
for  sensors  with  a  larger  bias  drift.  Furthermore,  even  if  this  debiasing  were  perfect,  it  still  does  not  address 
the  issue  of  residual  biases,  which  is  the  difference  between  the  true  sensor  bias  and  the  estimated  bias 
which  was  used  to  debias  the  sensors.  Residual  biases  must  be  treated  by  a  consider  analysis  or  by  the 
Schmidt-Kalman  filter  [2,  25],  Even  the  new  multi-billion  dollar  sensors  used  in  air  and  ground  tracking 
(e.g.,  SBX,  FBX,  TPY-2)  are  subject  to  these  issues,  so  the  older  sensors  that  make  up  a  large  portion  of  the 
space  surveillance  network  will  necessarily  exhibit  them  as  well. 

In  the  second  case,  the  association  is  unknown  (and  truth  is  unavailable)  and  both  the  association  and 
bias  must  be  solved  for  jointly.  This  is  a  much  more  difficult  problem,  and  it  is  generally  treated  with  what 
is  known  as  “pattern  matching”  [26-28].  Unlike  the  first  case,  which  relies  on  truth  or  known  associations, 
pattern  matching  can  be  performed  online  using  only  readily  available  data.  Because  of  this  the  joint  bias 
and  association  problem  is  the  key  bias  issue  needed  for  long  term  propagation,  and  the  remainder  of  this 
subsection  will  focus  on  it. 

The  problem  of  association  and  fusion  in  the  presence  of  biases  is  illustrated  in  Figure  7  (with  full  math¬ 
ematical  details  in  [26-28]).  Both  figures  show  estimates  of  object  positions  in  two-dimensional  space  with 
circles  representing  the  covariance  estimates.  The  large  covariances  correspond  to  previously-established 
tracks  that  have  been  predicted  into  the  held  of  view  of  the  next  sensor,  which  has  tracks  established  with  the 
smaller  covariances^.  Consider  first  the  situation  of  Figure  7(a),  where  bias  has  not  been  accounted  for.  Due 
to  the  overlap  of  the  “three-sigma”  error  circles,  four  tracks  from  each  sensor  will  be  incorrectly  assigned 
to  each  other,  while  three  tracks  from  each  sensor  remain  uncorrelated.  Compare  this  with  the  situation  of 

§  A  tracklet  is  a  short  track  segment  composed  from  about  10-12  sensor  measurements. 

^This  example  is  for  the  track-to-track  problem,  though  the  ideas  are  equally  applicable  for  measurement-to-track. 
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Figure  7:  (a)  The  effect  of  bias  on  data  association,  (b)  Data  association  after  bias  correction. 


Figure  7(b),  where  the  concurrent  bias  estimation  and  association  described  later  have  been  used  to  remove 
residual  sensor  biases  so  that  the  estimates  now  line  up  correctly.  In  this  case,  most  estimates  emanating 
from  common  objects  are  correctly  correlated,  with  two  estimates  that  do  not  correlate. 

The  bias  problem  for  handover,  track-to-track  association  ( correlation ),  and  fusion  described  above  has 
been  identified  as  a  major  obstacle  for  missile  defense  and  manifests  itself  even  for  the  newest  MDA  sensors; 
thus,  it  (along  with  the  analogous  measurement-to-track  problem)  will  certainty  present  difficulties  for  the 
UCT  problem  in  space  surveillance. 

We  conclude  this  subsection  with  a  brief  description  of  the  proposals  for  resolving  the  problem  of  joint 
association  and  bias  estimation.  The  problem  addressed  in  this  work  is  that  of  associating  two  classes  of 
objects  in  the  presence  of  a  bias  which  is  modeled  as  a  displacement  between  the  two  classes.  Here  it  is 
assumed  that  each  object  can  be  assigned  at  most  once  (this  can  be  generalized).  The  first  collection  of 
objects  is  enumerated  by  i  =  1, . . . ,  m,  and  the  second  by  j  =  l,...,n.  The  formulation  of  this  association 
problem  is  that  of  the  two-dimensional  inequality-constrained  assignment  problem,  but  with  an  unknown 
displacement  between  the  objects.  Inequality  constraints  are  used  since  not  all  objects  need  be  assigned. 

m  n 

Minimize^)  dT R~1d  +  ^  ^  Cij{d)xij 

i= 1  j= 1 

n 

Subject  To:  1  =  (16) 

3= 1 

m 

<1  (j  = 

i=l 

Xij  G  {0, 1}, 

where  d  =  (dx,  dy)  is  the  vector  of  prior  estimates  of  the  sensor  biases  and  R  is  the  correlation  matrix  for 
these  estimates.  The  correlation  costs  take  the  form  [21,  28] 

l 

(‘i.j ( d)  —  -(xj  -|-  dx  ( Ij'j  T  dy ) )  {Pa  T  Qjj  Sij  Sj ; )  {xi  T  dx  {yj  -f-  dy ) ) 

1 

T  —  In  ( det ( Pa  T  Qjj  Sij  'Yiji 

where  Pa  and  QVJ  are  the  covariances  of  the  track  state  estimates  xr  and  yj,  respectively,  while  St]  denotes 
the  cross-correlation  between  x,  and  y3 ,  which  may  arise  from  common  process  noise  or  common  priors  in 
the  sensor  estimation  filters. 

As  posed,  the  problem  (16)  falls  within  a  class  of  nonconvex  mixed-integer  nonlinear  programming 
problems  (MINLP)  [29-32].  The  techniques  for  solving  convex  MINLP  problems  include  outer  approxi- 
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mation,  generalized  Bender’s  decomposition,  extended  cutting  plane,  branch  and  bound,  and  branch  and  cut 
methods.  For  the  nonconvex  case,  one  attempts  to  extend  these  methods  by  using  convex  envelopes  or  by 
developing  convex  under  estimators  of  the  objective  function  and  constraints.  In  addition,  sampling  meth¬ 
ods  such  as  clustering  methods,  evolutionary  algorithms,  simulated  annealing,  and  tabu  search  may  also  be 
utilized. 

For  branch  and  bound,  the  lower  bound  that  one  must  compute  at  each  node  in  the  branch  and  bound  tree 
must  be  a  global  lower  bound.  To  achieve  this,  one  often  relaxes  the  discrete  variables  to  continuous  ones  and 
develops  convex  underestimators  for  the  objective  functions  and  constraints  that  facilitate  the  computation 
of  this  global  lower  bound.  An  example  of  this  is  the  a-branch  and  bound  method  of  Adjiman,  Androulakis, 
and  Floudas  [33],  which  can  be  applied  to  the  current  problem. 

The  branch  and  bound  framework  developed  herein  does  not  follow  these  approaches  for  MINLP  prob¬ 
lems.  The  constraints  arc  already  convex  in  that  they  arc  affine.  The  objective  function  is  nonconvex; 
however,  we  can  develop  a  convex  (affine)  lower  approximation  to  the  objective  function  by  using 

bij  =  ^  ln(det(Pjj  +  Qn  -  Sij  -  Sji ))  -  7^  for  i  /  0  and  j  f  0,  and  (17) 

bij  =  Cij  for  i  =  0  or  j  =  0, 

both  of  which  are  independent  of  the  displacement  d.  Furthermore,  bij  <  c%3  (df  thereby  providing  a  basis 
for  the  global  lower  bound  needed  in  the  branch  and  bound  procedure.  We  can  make  effective  use  of  an 
assignment  solver  so  there  is  no  need  for  relaxation  of  the  discrete  variables. 

For  space  surveillance,  a  multidimensional  (N  >  2)  version  of  the  joint  bias  estimation  and  correlation 
problem  must  be  formulated  and  likelihood  ratios  derived.  Such  algorithms  would  be  based  on  Lagrangian 
relaxation  and  an  anytime  branch  and  bound  algorithm  similar  to  that  developed  for  the  two-dimensional 
assignment  problem  with  appropriate  upper  and  lower  bounds  on  optimality. 

1. 1.2.4  Uncertainty  in  Correlation  The  goal  of  this  algorithm  is  the  development  of  an  approximation 
to  the  probability  of  correct  association.  This  has  been  achieved  for  the  sequential  k-best  approximation  to 
MHT  based  either  on  a  k-best  approach  [34]  or  Markov  Chain  Monte  Carlo  (MCMC)  methods  [35],  but  not 
yet  for  the  full  multidimensional  assignment  problem. 

In  the  k-best  approach,  an  MHT  system  directly  maintains  a  set  77 /.  of  k  “complete”  data  association  hy¬ 
potheses  Hi  £  Hk,  i  =  1 .... .  k  within  a  sliding  window  of  sensor  data  frames  to  approximate  the  optimal 
solution  that  maximizes  the  data  association  problem.  Each  alternative  solution  represents  a  different  com¬ 
plete  data  association  hypothesis  Hi.  Let  Hk  =  {Hi,  H2,  ■  .  ■ ,  Hk}  denote  a  set  of  k  ranked  solutions,  such 
that  H\  corresponds  to  the  “best”  solution  returned  from  the  assignment  solver,  H>  to  the  “second-best”, 
etc.;  that  is, 

chAn)  <ch2{N)  <...<  cHk(N), 

where  the  costs  are  defined  as  negative  log-likelihood  terms.  As  long  as  k  is  sufficiently  large,  it  is  the  case 
that 

k 

^>t{Hq\Zn)  +  Fi’(Hi\ZN)  «  1, 

i= 1 

and  therefore 

Fr(H  \ZN)  =  PlWZ^  ~  Pr mZN) 

1  Pi\H0\ZN)  +  E  •=!  Pr^-I-Z*) ' 
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Using  (13)  and  (14),  we  obtain  [21,  36] 


Pr(fli) 


*-chAN) 


l  +  £j=1e 


which  gives  the  probability  of  a  hypothesis  in  terms  of  its  cost  CHi- 

Now,  let  Xij  denote  the  probabilistic  data  association  event  that  measurement  associates  with  y3. 
Then,  the  probability  Pr (xij)  °f  this  event  is  given  by 


Pr(Xii)«  E  Pr(^E 

HeHn 


where  7itJ  C  74/,.  denotes  the  subset  of  those  data  association  hypotheses  from  the  set  74/,  that  postulate  the 
event  x%j- 


1.1.2.5  Results  Though  there  is  still  much  work  to  be  done  in  developing  an  MHT  for  space  surveillance, 
some  initial  investigations  have  been  conducted  under  this  effort.  Here  we  will  present  results  from  two 
breakup  scenarios,  which  are  described  in  Table  2.  For  each  scenario  we  will  compare  MHT  results  using 
different  nonlinear  filters  including  (a)  a  single  Gaussian  (UKF),  (b)  a  course  Gaussian  sum,  (c)  a  fine 
Gaussian  sum,  and  (d)  the  non-statistical  method  employed  in  the  AFSPC  Astrodynamic  Standards  ROTAS 
program.  The  number  of  terms  in  the  Gaussian  sum  is  given  by  TV,  with  TV  =  1  indicating  a  single  Gaussian, 
and  TV  =  it  indicating  the  non-statistical  (ROTAS)  method. 


Table  2:  Description  of  two  breakup  scenarios. 


Scenario  One 

Scenario  Two 

Regime 

LEO 

LEO 

Bias,  drag,  maneuvers 

Not  taken  into  account 

No.  of  tracklets 

27 

42 

Duration 

6.5  hours 

24  hours 

No.  of  sensors 

4 

1 

Figure  8  compares  MHT  results  between  pairs  of  methods  for  the  two  scenarios.  The  axes  arc  indices 
into  the  tracklets  and  each  pixel  tells  whether  the  corresponding  pair  of  tracklets  arc  in  the  same  track  in 
both  methods  (black),  neither  method  (white),  method  one  only  (red),  or  method  two  only  (green).  Note 
that  black  and  white  pixels  indicate  agreement  between  the  methods,  while  red  and  green  pixels  indicate 
disagreement.  Table  4  supplements  the  information  in  the  figures  with  the  metrics  defined  in  Table  3. 


Table  3:  Definitions  of  the  track  and  summary  metrics. 


Track  Metric  T) 
Track  Metric  T2 
Track  Metric  T& 

Summary  Metric  S 


Number  of  tracks  in  the  first  method 
Number  of  tracks  in  the  second  method 

Number  of  tracks  that  appear  in  both  methods  with  the  same  as¬ 
sociated  tracklets 

Percentage  of  tracklet  pairs  for  which  both  methods  agree  on 
whether  or  not  they  emanate  from  a  common  object 


The  left  three  plots  in  Figure  8  and  the  left  half  of  Table  4  show  the  results  for  Scenario  One.  For  this 
scenario,  the  single  Gaussian  (UKF),  coarse  Gaussian  sum  (TV  =  17),  and  fine  Gaussian  sum  (TV  =  56) 
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Scenario  One 


Scenario  Two 


N  =  it  vs.  N  =  56  N  =  it  vs.  N  =  121 


Figure  8:  Plots  comparing  MHT  results  using  different  nonlinear  filters. 
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Table  4:  Track  and  summary  metrics  (defined  in  Table  3)  for  Scenario  One  (left)  and  Scenario  Two  (right) 


N\ 

n2 

C Ti,T2,Tb ) 

5 

tv, 

n2 

C Ti,T2,Tb ) 

5 

1 

56 

(18,18,18) 

1.000000 

1 

121 

(26,24,11) 

0.973923 

17 

56 

(18,18,18) 

1.000000 

56 

121 

(24,24,22) 

0.995465 

it 

56 

(18,18,13) 

0.980796 

it 

121 

(37,24,6) 

0.976190 

all  give  the  same  answers.  Specifically,  all  three  methods  produce  the  exact  same  18  tracks.  From  this 
we  can  conclude  that  in  breakup  scenarios  with  short  time  intervals  between  updates,  propagating  a  state 
and  covariance  using  the  UKF  is  robust  and  sufficient  -  the  higher  order  representations  have  validated  this 
conclusion.  In  contrast,  the  non-statistical  approach  also  produces  18  tracks,  but  only  13  of  those  tracks  arc 
the  same  as  the  tracks  produced  with  the  statistical  methods.  This  indicates  that  even  over  a  short  scenario 
the  non-statistical  approach  is  insufficient  and  need  not  produce  the  correct  solution. 

The  right  three  plots  in  Figure  8  and  the  right  half  of  Table  4  show  the  results  for  Scenario  Two.  Here 
we  see  a  reasonable  amount  of  disagreement  between  the  single  Gaussian  (UKF)  and  fine  Gaussian  sum 
(N  =  121),  but  almost  exact  agreement  between  the  coarse  Gaussian  sum  (N  =  56)  and  the  fine  Gaussian 
sum.  This  “convergence”  as  the  fidelity  of  the  filter  increases  provides  evidence  that  advanced  filtering 
techniques  such  as  adaptive  Gaussian  sum  filters  can  be  used  to  produce  the  correct  solutions.  It  is  interesting 
to  look  a  bit  closer  at  the  form  of  the  discrepancy  between  the  two  Gaussian  sum  filters.  The  coarse  Gaussian 
sum  associates  tracklets  (17,  36)  and  (23,  35),  while  the  fine  Gaussian  sum  switches  these  pairs  to  (17,  35) 
and  (23,  36)  -  the  only  difference  in  results  between  these  two  methods  is  a  single  association  decision. 
Moving  to  the  non-statistical  approach,  we  see  a  stark  contrast.  The  non-statistical  approach  produced  37 
tracks,  indicating  that  it  had  trouble  associating  objects  together  (we  can  also  see  this  by  the  large  number  of 
green  pixels  in  the  figure).  In  particular,  it  made  very  few  associations  between  tracklets  that  were  separated 
by  long  time  intervals  -  the  statistical  methods  did  not  have  this  problem. 

Finally,  we  mention  that  the  results  presented  here  only  focus  on  the  best  hypothesis  from  each  method. 
In  cases  where  the  best  hypothesis  has  less  than  100%  of  the  total  likelihood,  the  ambiguity  must  also  be 
considered  in  order  to  get  “the  full  picture.”  Even  in  these  relatively  small  scenarios  there  was  significant 
ambiguity,  though  for  simplicity  the  ambiguity  is  not  included  in  the  presented  results. 

1. 1.2.6  Status  The  small  breakup  scenarios  in  LEO  discussed  above  were  presented  at  the  AFOSR  re¬ 
view  meeting  on  Sept.  11,  2010  and  at  the  DARPA  meeting  on  UCTs  at  Washington  DC  on  Nov.  17,  2010. 
Extensions  to  GEO  and  the  joint  association  and  bias  estimation  problem  are  ongoing. 

1.1.3  Unresolved  Closely  Spaced  Objects 

Much  of  this  section  discusses  algorithm  components  that  have  been  developed  under  a  Phase  II  SBIR, 
“Hierarchical  Image  Processing  for  Closely  Spaced  Object  (CSO)  Resolution”.  The  core  research  problem 
addressed  in  this  paid  of  the  project  is  that  of  augmenting  measurement  generation  algorithms  to  supply  high- 
quality  measurement  covariances  and  estimates  of  the  uncertainty  in  object  count  in  a  small  region  of  space, 
called  resolution  uncertainty,  for  use  in  generating  accurate  orbit  state  covariances  and  data  association 
ambiguity  information. 

The  innovative  paid  of  our  approach  is  the  use  of  multiframe  ideas  in  two  manners.  First,  we  propose 
to  use  multiple  images  from  the  sensor  to  generate  high-quality  individual  images  using  image  superreso¬ 
lution  techniques.  Second,  given  a  set  of  images,  either  directly  from  the  sensor  or  preprocessed  using  the 
above  techniques,  we  will  use  multiple-frame,  multiple-hypothesis  expectation  maximization  (EM)-based 
measurement  generation  to  solve  the  problem  of  targets  appealing  in  some  frames  while  not  appealing  in 
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others,  the  so-called  stochastic  resolution  problem.  A  combination  of  such  methods  promises  to  substan¬ 
tially  improve  the  resolution  and  tracking  of  closely  spaced  objects. 

1. 1.3.1  Background  on  Image  Processing  The  field  of  image  processing  for  astronomy  and  space 
surveillance  has  involved  a  significant  amount  of  research  focused  on  the  problem  of  improving  image 
quality  prior  to  extracting  target  measurement  data.  Single-frame  image  restoration  methods  involve  various 
techniques  intended  to  reduce  the  noise  in  the  image,  reduce  the  blurring  that  results  from  optical  diffraction 
and  atmospheric  turbulence,  or  both.  In  general,  noise  reduction  and  blur  compensation  are  at  odds;  most 
methods  that  reduce  blur  will  tend  to  amplify  noise,  and  vice  versa.  These  techniques  include  basic  spa¬ 
tial  filtering  [37],  deconvolution  approaches  [38-40],  variations  on  the  CLEAN  algorithm  [41],  and  newer 
techniques  based  on  wavelet  transform  filtering  [42]. 

All  of  the  above  methods  for  image  restoration  are  intended  to  function  using  a  single  input  frame  of  data 
and,  usually,  some  form  of  prior  knowledge  (such  as  the  image  noise  statistics  and  the  point  spread  function) 
to  generate  a  restored  scene.  One  of  the  core  innovations  in  this  work  is  to  broaden  the  view  of  the  problem 
to  include  multiframe  methods,  in  which  data  from  more  than  one  collected  image  may  be  used  to  generate  a 
result.  This  additional  data  permits  both  resolution  enhancement  and  consideration  of  measurement  stability 
in  the  presence  of  unresolved  closely  spaced  objects.  Multiframe  image  restoration  methods  in  the  literature 
range  from  shift-and-add  formulations  with  anti-aliasing  through  hybrid  approaches  such  as  the  Drizzle 
algorithm  [43]  all  the  way  up  through  to  image  superresolution  techniques. 

Most  of  the  image  restoration  algorithms  discussed  in  the  literature  arc  pure  image  transforms;  that  is, 
they  take  as  input  one  or  more  images  and  produce  a  higher-quality  output  image.  This  alone  is  insufficient 
for  the  puiposes  of  tracking;  we  must  convert  the  raw  image  data  into  a  sequence  of  measurements  that  can 
be  correlated  and  fused  to  produce  orbit  state  estimates.  This  process  can  be  challenging  during  events  such 
as  on-orbit  collisions,  where  potentially  large  numbers  of  objects  may  reside  within  the  field  of  view  of  a 
sensor.  Prior  work  by  Numerica  on  tracking  from  infrared  sensor  data  [44]  has  shown  that  multiple-frame, 
multiple-hypothesis  expectation  maximization  (EM)-based  measurement  generation  has  the  ability  to  solve 
the  stochastic  resolution  problem  by  considering  information  from  several  consecutive  frames  of  data  to 
determine  the  correct  measurement  count  in  the  scene.  Such  algorithms  stabilize  the  observed  object  scene 
in  the  presence  of  potentially  significant  numbers  of  closely  spaced  objects  and  arc  discussed  below. 

1.1.3.2  Proposed  Implementation  Consider  the  case  of  an  imaging  sensor  viewing  one  or  more  objects 
as  point  sources  (that  is,  the  sensor  resolution  is  insufficient  to  observe  the  objects  in  more  than  one  pixel). 
The  observed  image  Y  is  a  function  of  object  position  and  magnitude,  plus  noise: 

K 

Y  =  ]Tg^)  +  E, 

1=1 

where  0r  is  a  vector  [//„,; ,  pvi ,  I,}'  containing  the  centroid  focal  plane  coordinates  u.  v  and  intensity  magni¬ 
tude  I  for  object  i.  Further,  E  is  the  additive  noise  present  in  the  image,  and  G,  is  the  (nonlinear)  imaging 
function.  We  would  like  to  find  an  estimate  Oi  for  each  object  i  that  maximizes  the  likelihood  of  the  ob¬ 
served  image.  It  may  be  observed  that  the  structure  of  this  problem  is  similar  to  that  of  the  multiframe 
image  restoration  problem  presented  by  Elad  and  Feuer  [45];  this  analogy  is  intentional.  We  will  discuss  the 
combined  image  restoration  and  measurement  extraction  problem  below. 

In  the  special  case  where  the  imaging  function  can  be  approximated  by  a  Gaussian  blur  function  resulting 
from  sensor  limitations  and  atmospheric  turbulence,  we  can  be  more  precise  in  the  formulation  of  the  above 
equation.  We  collapse  the  2D  image  data  into  a  column  vector  for  ease  of  manipulation.  For  a  single 
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measured  pixel  yj,  the  relationship  between  the  point  sources  and  the  pixel  intensity  is  given  by 

1  1  (  |  (V-Vyi)2  2p(u-/xui)(u-p„i)  . 

=e  2  '-p'2  v  °2  °2  <TuCTi'  dudv  +  E.i 

9  J 

P2 

for  an  arbitrary  Gaussian  blur  function  B.  If,  as  is  typical,  the  blur  function  is  spherical  (that  is,  B  =  n2I), 
then  p  =  0  in  the  above  equation.  This  form  is  separable,  and  therefore  yj  may  be  expressed  as  a  closed-form 
solution  in  terms  of  the  error  function  for  computational  efficiency. 

We  now  have  a  model  for  the  intensity  of  each  pixel  in  the  image  that  we  can  use  to  formulate  a  nonlinear 
estimation  problem  for  the  point  source  locations  and  intensities.  This  nonlinear  estimation  problem,  how¬ 
ever,  tends  to  be  very  poorly  conditioned  and  have  local  minima.  Good  starting  guesses  for  the  point  source 
locations  arc  required  if  the  estimator  is  to  converge  to  the  correct  solution.  Also,  it  is  clear  that  this  model 
requires  some  prior  estimate  of  the  number  K  of  point  sources  to  be  sought  in  the  image.  In  some  cases  this 
estimate  will  be  available  from  the  sensor  tasking  (if  the  objects  in  the  field  of  view  are  well-characterized); 
however,  this  is  not  typically  the  case  for  objects  resulting  from  new  space  events.  The  key  to  solving  this 
problem  lies  in  the  multiple  hypothesis  approach.  We  have  constructed  a  range  of  alternate  approaches  to 
generating  high-probability  prior  object  counts  and  locations.  Then,  given  these  priors,  we  can  construct 
several  different  measurement  hypotheses  for  the  image. 

We  have  observed  that  the  most  effective  method  for  generating  a  measurement  prior,  if  the  input  datarate 
supports  use  of  the  technique,  is  that  of  first  constructing  a  synthetic  high-resolution  image  from  a  sequence 
of  low-resolution  images  (multiple  frame  image  restoration).  In  this  case,  the  formulation  of  the  problem  is 
now 


K 


Vi  =  *52 


i=  1 


2 Troy  (7 v  \/l 


Y  =  H(X)  +  E, 

I< 

x=x:  ga 

i=i 

with  a  lineal-  transform  H  representing  blur,  warp,  and  decimation  functions  between  a  theoretical  high- 
resolution  image  X  and  the  measured  low-resolution  images.  We  accumulate  sufficient  measured  images 
to  solve  the  sparse  linear  estimation  problem  for  the  high-resolution  image,  and  then  use  standard  peak 
extraction  techniques  to  obtain  a  prior  for  the  nonlinear  parameter  estimation  problem  described  above. 

While  the  use  of  multiframe  image  restoration  for  astronomy  applications  is  not  novel,  our  unique  ap¬ 
proach  combines  the  image  restoration  problem  with  source  position  estimation  for  tracking  in  the  presence 
of  unresolved  closely  spaced  objects.  Figure  9  shows  the  multiple  frame  image  restoration  used  in  conjunc¬ 
tion  with  model-based  measurement  generation  to  achieve  maximum  point  source  resolving  power.  Panel  (a) 
shows  a  raw  measured  image  (single -band  data  in  false  color)  of  a  set  of  point  sources  at  locations  marked 
with  white  ‘+’  symbols,  along  with  measurements  and  covariances  (blue  V  and  ellipses)  generated  from  a 
simple  cluster  and  centroid  approach.  Panel  (b)  focuses  on  a  small  subset  of  the  image  containing  several 
unresolved  point  sources.  Panel  (c)  is  a  synthetic  high-resolution  image  of  the  same  scene  generated  from 
a  sequence  of  images  such  as  (a),  with  measurements  constructed  from  the  nonlinear  optimization  problem 
described  above.  The  advanced  measurement  generation  algorithms  have  the  ability  to  resolve  point  sources 
with  much  greater  efficacy  than  would  be  possible  from  a  single-frame  approach. 

Using  multiple  consecutive  frames  of  data,  we  consider  whether  the  evolution  of  measurements  over 
time  produces  a  consistent  track  picture.  In  order  to  solve  this  problem  within  an  orbit  determination  context, 
we  require  two  additional  items  of  information:  the  covariance  associated  with  each  of  the  estimated  object 
locations,  and  the  uncertainty  of  the  object  resolution  hypothesis.  The  measurement  covariance  is  a  function 
of  the  signal-to-noise  ratio  for  the  object  being  observed.  If  there  are  multiple  objects  within  a  confusable 
distance  (overlapping  blur  regions)  of  one  another,  then  the  error  associated  with  the  centroid  estimate  is 
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also  a  function  of  the  distance  to,  and  brightness  of,  the  nearby  objects.  Failure  to  account  for  these  effects 
can  result  in  an  incorrect  measurement  covariance  estimate.  We  present  here  the  development  for  object 
resolution  uncertainty  given  two  objects  within  a  confusion  distance  from  one  another;  development  of  the 
more  general  form  is  anticipated  as  part  of  this  research  program.  Suppose  we  have: 

Y  =  G(0)+E~iV(G(0),r72l); 


that  is,  the  intensity  noise  in  the  image  is  Gaussian  and  characterized  by  noise  power  o\. 

Let  us  consider  the  hypotheses  that  6q  =  0  (the  data  consists  of  noise  only),  0\  =  [60bj\ ;  0]  (the  data 
consists  of  object  1  only),  O2  =  [0;  0obj2]  (the  data  consists  of  object  2  only),  and  #3  =  [Q0bj  1 ;  0obj2]  (both 
objects  must  be  present  to  explain  the  data).  The  generalized  likelihood  of  each  of  these  hypotheses  can  then 
be  expressed  as 


fi{y)  =  (2ir)-n/2a~n  exp 


These  likelihoods  allow  us  to  resolve  the  hypotheses  described  above.  Critical  for  the  resolution  question 
is  the  ratio  of  the  likelihood  that  the  data  is  explained  by  one  object  versus  two.  If  the  algorithm  has 
hypothesized  the  presence  of  two  objects,  but  it  is  more  likely  that  the  data  can  be  explained  by  only  one 
object,  then  we  either  (a)  really  do  only  have  one  object,  or  (b)  have  two  objects  that  cannot  be  distinctly 
resolved.  These  likelihood  ratios  then  form  the  basis  for  incorporating  measurement  hypothesis  probabilities 
into  the  tracking  problem. 


1. 1.3.3  Conclusion  We  have  discussed  two  innovative  approaches  that  combine  to  address  tracking 
problems  for  closely  spaced  objects,  both  using  multiframe  ideas,  but  in  two  subtly  different  ways.  First, 
we  discussed  using  multiple  images  from  the  sensor  to  generate  high-quality  individual  images  using  image 
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superresolution  techniques.  Second,  given  a  set  of  images,  either  directly  from  the  sensor  or  preprocessed 
using  the  above  techniques,  we  showed  how  to  use  multiple -frame,  multiple-hypothesis  expectation  max¬ 
imization  (EM)-based  measurement  generation  to  solve  the  problem  of  targets  appealing  in  some  frames 
while  not  appealing  in  others,  giving  rise  to  a  possible  solution  to  both  stochastic  resolution  and  resolution 
uncertainty  problems.  We  have  constructed  a  prototype  of  the  combined  multiframe  approach;  future  work 
will  focus  on  use  of  this  formulation  to  construct  measurement  covariances  and  resolution  uncertainties  for 
tracking. 

1.1.3.4  Status  This  work  is  ongoing.  The  processing  of  data  in  GEO  has  not  yet  started. 

1.2  Accomplishments  /  New  Findings 

The  primary  objective  of  this  program  is  to  perform  the  necessary  basic  research  to  support  the  development 
of  a  statistical,  multiple  hypothesis  tracker  (MHT)  for  space  surveillance.  Such  a  MHT  framework  can  serve 
as  the  next  generation  space  surveillance  system  to  maintain  the  space  catalog,  to  identify  uncorrelated 
tracks,  and  to  support  conjunction  analysis  and  sensor  resource  management.  Key  components  in  such 
a  system  include  a  consistent  characterization  of  uncertainty,  physical  modeling,  multiple  model  filtering, 
and  the  association  problem  of  determining  which  tracklets/measurements  emanate  from  which  object.  To 
achieve  a  consistent  characterization  of  uncertainty,  Numerica  has  developed  an  adaptive  Gaussian  sum 
filter  which  correctly  represents  and  propagates  uncertainties  and  adaptively  selects  the  correct  the  number 
of  Gaussians  in  the  mixture.  Realtime  online  metrics  support  the  coarsening  and  refining  of  the  filter  to 
maintain  consistent  uncertainty.  A  sliding  window  batch  estimation  filter  has  also  been  developed  and  shown 
to  provide  an  accurate  evaluation  of  the  prediction  error  critical  for  collect  anomaly  detection  and  resolution 
of  UCTs.  Several  papers  have  been  published  on  these  uncertainty  management  and  propagation.  Numerica 
has  also  acquired  orbital  data  from  JFCC  SPACE  and  processing  of  this  data  has  commenced  using  a  new 
prototype  MHT  specifically  suited  to  space  surveillance  in  conjunction  with  the  aforementioned  techniques 
on  uncertainty  management. 

2  Personnel  Supported 

a.  PI:  Aubrey  B.  Poore 

b.  Colleagues  at  Numerica:  Joshua  T.  Horwood,  Nathan  D.  Aragon,  and  Scott  Danford 


3  Publications 

Journal  papers  submitted  for  publication  under  this  effort  are  cited  as  References  [5,  17]  in  the  bibliography. 
Conference  proceedings  papers  produced  under  this  effort  are  cited  as  References  [15,  16,  18,  19]  in  the 
bibliography. 

4  Participation  /  Presentations  at  Meetings,  Conferences,  Seminars,  etc. 

1.  SPIE  Defense,  Security,  and  Sensing  2010:  Signal  and  Data  Processing  of  Small  Targets.  Orlando, 
FL.  Apr.  5-9,  2010.  Presented  paper  [15]. 

2.  George  H.  Born  Symposium.  Boulder,  CO.  May  13-14,  2010. 

3.  Kyle  T.  Alfriend  Astrodynamics  Symposium.  Monterey,  CA.  May  17-19,  2010.  Presented  paper  [18]. 
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4.  Advanced  Maui  Optical  and  Space  Surveillance  Technologies  Conference  2010.  Wailea,  HI.  Sept.  15- 
17,  2010.  Presented  paper  [19]. 

5.  DARPA  UCT  meeting.  Washington,  DC.  Nov.  17,  2010. 

6.  21st  AAS/AIAA  Space  Flight  Mechanics  Meeting.  New  Orleans,  LA.  Feb.  14-17,  2011.  Presented 
paper  [16]. 


5  New  Discoveries,  Inventions,  or  Patents 

No  patents  resulted  from  this  effort. 
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